1. Background

The HVT package is a collection of R functions to facilitate building topology preserving maps for rich multivariate data analysis. Tending towards a big data preponderance, a large number of rows. A collection of R functions for this typical workflow is organized below:

1. Data Compression: Vector Quantization (VQ), HVQ (Hierarchical Vector Quantization) using means or medians. This step compresses the rows (long data frame) using a compression objective.

2. Data Projection: Dimension projection of the compressed cells to 1D, 2D and Interactive surface plot with the Sammons Non-linear Algorithm. This step creates a topology preserving map (also called as mathematical embeddings) coordinates into the desired output dimension.

3. Tessellation: Create cells required for object visualization using the Voronoi Tessellation method, package includes heatmap plots for Hierarchical Voronoi Tessellations (HVT). This step enables data insights, visualization, and interaction with the topology preserving map. Useful for semi-supervised tasks.

4. Scoring: Scoring data sets and recording their assignment using the map objects from the above steps, in a sequence of maps if required.

5. Temporal Analysis and Visualization: A Collection of functions that leverages the capacity of the HVT package by analyzing time series data for its underlying patterns, calculation of transitioning probabilities and the visualizations for the flow of data over time.

What’s New?

• The new update focuses on the integration of time series capabilities into the HVT package by extending its foundational operations to time series data which is emphasized in this vignette.

• The new functionalities are introduced to analyze underlying patterns and trends within the data, providing insights into its evolution over time and also offering the capability to analyze the movement of the data by calculating its transitioning probability and creates elegant plots and GIFs.

Below are the new functions and its brief descriptions:

• plotStateTransition: Provides the time series flowmap plot.
• getTransitionProbability: Provides a list of transition probabilities.
• reconcileTransitionProbability: Provides plots and tables for comparing transition probabilities calculated manually and from markovchain function.
• plotAnimatedFlowmap: Creates flowmaps and animations for both self state and without self state scenarios.

2. Experimental setup

The Lorenz attractor is a three-dimensional figure that is generated by a set of differential equations that model a simple chaotic dynamic system of convective flow. Lorenz Attractor arises from a simplified set of equations that describe the behavior of a system involving three variables. These variables represent the state of the system at any given time and are typically denoted by (x, y, z). The equations are as follows:

\[ dx/dt = σ*(y-x) \] \[ dy/dt = x*(r -z)-y \] \[ dz/dt = x*y-β*z \] where dx/dt, dy/dt, and dz/dt represent the rates of change of x, y, and z respectively over time (t). σ, r and β are constant parameters of the system, with σ(σ = 10) controlling the rate of convection, r(r=28) controlling the difference in temperature between the convective and stable regions, and β(β = 8/3) representing the ratio of the width to the height of the convective layer. When these equations are plotted in three-dimensional space, they produce a chaotic trajectory that never repeats. The Lorenz attractor exhibits sensitive dependence on initial conditions, meaning even small differences in the initial conditions can lead to drastically different trajectories over time. This sensitivity to initial conditions is a defining characteristic of chaotic systems.

In this notebook, we will use the Lorenz Attractor Dataset. This dataset contains 200,000 (Two hundred thousand) observations and 5 columns. The dataset can be downloaded from here.

The dataset includes the following columns:

• X: The X-coordinate of a point in the Lorenz attractor. It represents the measure of strength of the circulation of fluid in the convective flow
• Y: The Y-coordinate of a point in the Lorenz attractor. It represents the measure of the temperature of the fluid at a given point in space, but it is the horizontal temperature distribution
• Z: The Z-coordinate of a point in the Lorenz attractor. It represents the measure of the temperature of the fluid at a given point in space, but it is the vertical temperature distribution
• U: U represents the velocity or speed of the system at a particular point in the attractor
• t: The time variable (discrete time steps which are incremented by approximately 0.0002500012500062 units of time) associated with each point in the Lorenz attractor. It indicates how much time has elapsed since the beginning of the simulation. Each value represents a specific point in time, and they are spaced apart by a fixed time interval.

3. Notebook Requirements

This chunk verifies the installation of all the necessary packages to successfully run this vignette, if not, installs them and attach all the packages in the session environment.

list.of.packages <- c("dplyr", "kableExtra", "plotly", "purrr", "data.table", "gridExtra", "grid", "reactable", "reshape", "tidyr", 
                      "stringr", "DT", "knitr", "feather")

new.packages <-
  list.of.packages[!(list.of.packages %in% installed.packages()[, "Package"])]
if (length(new.packages))
  install.packages(new.packages, dependencies = TRUE, verbose = FALSE, repos='https://cloud.r-project.org/')
invisible(lapply(list.of.packages, library, character.only = TRUE))

# Sourcing required code scripts for HVT
script_dir <- "../R"
r_files <- list.files(script_dir, pattern = "\\.R$", full.names = TRUE)
invisible(lapply(r_files, function(file) { source(file, echo = FALSE); }))

4. Data Understanding

Here, we load the data. Let’s explore the Lorenz Attractor Dataset. For the sake of brevity, we are displaying only the first ten rows.

file_path <- ("./sample_dataset/lorenz_attractor.feather")
dataset <- read_feather(file_path) %>% as.data.frame()
dataset <- dataset %>% select(X,Y,Z,U,t)
dataset$t <- round(dataset$t, 5)
displayTable(head(dataset, 10))

Now, let’s try to visualize the Lorenz attractor (overlapping spirals) in 3D Space.

data_3d <- dataset[sample(1:nrow(dataset), 1000), ]
plot_ly(data_3d, x= ~X, y= ~Y, z = ~Z) %>% add_markers( marker = list(
                          size = 2,
                          symbol = "circle",
                          color = ~Z,
                          colorscale = "Bluered",
                          colorbar = (list(title = 'Z'))))

Now, let’s have a look at the structure of the Lorenz Attractor dataset.

str(dataset)

## 'data.frame':    200000 obs. of  5 variables:
##  $ X: num  0 0.0025 0.00499 0.00747 0.00995 ...
##  $ Y: num  1 1 1 0.999 0.999 ...
##  $ Z: num  20 20 20 20 19.9 ...
##  $ U: num  0 0.0005 0.001 0.0015 0.002 ...
##  $ t: num  0 0.00025 0.0005 0.00075 0.001 0.00125 0.0015 0.00175 0.002 0.00225 ...

Data distribution This section displays four objects.

Variable Histograms: The histogram distribution of all the features in the dataset.

Box Plots: Box plots for all the features in the dataset. These plots will display the median and Interquartile range of each column at a panel level.

Correlation Matrix: This calculates the Pearson correlation which is a bivariate correlation value measuring the linear correlation between two numeric columns. The output plot is shown as a matrix.

Summary EDA: The table provides descriptive statistics for all the features in the dataset.

• variable: The features/columns of the dataset
• min: Minimum value of that feature/column
• 1st Quartile: The value that splits the lower 25% of the data when arranged in ascending order
• median: Middle value in the ascendingly ordered dataset
• mean: Sum of all values in the dataset divided by the total number of values
• sd: Measure of the dispersion of dataset relative to its mean.
• 3rd Quartile: The value that splits the lower 75% of the data when arranged in ascending order
• max: Maximum value of that feature/column
• hist: The basic barchart of the data distribution of a feature/column
• n_row: Number of rows for that feature/column
• n_missing: Number of missing values/NAs for that feature/column

Time Series Plots: Plots of all features (including time) against the time column.

It uses an inbuilt function called edaPlots to display the above-mentioned four objects.

NOTE: The input dataset should be a data frame object and the columns should be only numeric type.

edaPlots(dataset, time_column = "t", output_type = "timeseries", n_cols = 5)

edaPlots(dataset, output_type = 'summary', n_cols = 5)

edaPlots(dataset, output_type = 'histogram', n_cols = 5)

edaPlots(dataset, output_type = 'boxplot', n_cols = 5)

edaPlots(dataset, output_type = 'correlation', n_cols = 5)

Train - Test Split

Let us split the dataset into train and test. We will orderly select 80% of the data as train and the remaining as test.

noOfPoints <- dim(dataset)[1]
trainLength <- as.integer(noOfPoints * 0.8)
trainDataset <- dataset[1:trainLength,]
testDataset <- dataset[(trainLength+1):noOfPoints,]
rownames(testDataset) <- NULL

displayTable(head(trainDataset, 10))

str(trainDataset)

## 'data.frame':    160000 obs. of  5 variables:
##  $ X: num  0 0.0025 0.00499 0.00747 0.00995 ...
##  $ Y: num  1 1 1 0.999 0.999 ...
##  $ Z: num  20 20 20 20 19.9 ...
##  $ U: num  0 0.0005 0.001 0.0015 0.002 ...
##  $ t: num  0 0.00025 0.0005 0.00075 0.001 0.00125 0.0015 0.00175 0.002 0.00225 ...

edaPlots(trainDataset, time_column = "t", output_type = "timeseries", n_cols = 5)

edaPlots(trainDataset, output_type = 'summary', n_cols = 5)

edaPlots(trainDataset, output_type = 'histogram', n_cols = 5)

edaPlots(trainDataset, output_type = 'boxplot', n_cols = 5)

edaPlots(trainDataset, output_type = 'correlation', n_cols = 5)

displayTable(head(testDataset,10))

str(testDataset)

## 'data.frame':    40000 obs. of  5 variables:
##  $ X: num  16.1 16.1 16 16 16 ...
##  $ Y: num  13.7 13.6 13.6 13.5 13.5 ...
##  $ Z: num  39.6 39.6 39.7 39.7 39.7 ...
##  $ U: num  9.89 9.89 9.89 9.89 9.89 ...
##  $ t: num  40 40 40 40 40 ...

edaPlots(testDataset, time_column = "t", output_type = "timeseries", n_cols = 5)

edaPlots(testDataset, output_type = 'summary', n_cols = 5)

edaPlots(testDataset, output_type = 'histogram', n_cols = 5)

edaPlots(testDataset, output_type = 'boxplot', n_cols = 5)

edaPlots(testDataset, output_type = 'correlation', n_cols = 5)

5. Model Training and Visualization

We will use the trainHVT function to compress our dataset while preserving essential features.

Model Parameters

• Number of Cells at each Level = 100
• Maximum Depth = 1
• Quantization Error Threshold = 0.1
• Error Metric = Max
• Distance Metric = Manhattan / L1_Norm
• Dimension Reduction Metric = Sammon

NOTE: The compression takes place only for the X, Y and Z coordinates and not for U(velocity) and t(Timestamp). After training & Scoring, we merge back the U and t columns with the dataset.

hvt.results <- trainHVT(
  trainDataset[,-c(4:5)],
  n_cells = 100,
  depth = 1,
  quant.err = 0.1,
  normalize = TRUE,
  distance_metric = "L1_Norm",
  error_metric = "max",
  quant_method = "kmeans",
  dim_reduction_method = "sammon")

## Initial stress        : 0.00301
## stress after  10 iters: 0.00155, magic = 0.500
## stress after  20 iters: 0.00154, magic = 0.500

Let’s check out the compression summary.

displayTable(data = hvt.results[[3]]$compression_summary)

NOTE: Based on the provided table, it’s evident that the ‘percentOfCellsBelowQuantizationErrorThreshold’ value is zero, indicating that compression hasn’t taken place for the specified number of cells, which is 100. Typically, we would continue increasing the number of cells until at least 80% compression occurs. However, in this vignette demonstration, we’re not doing so, because the plots generated from temporal analysis functions would become cluttered and complex, making explanations less clear.

Let’s check out the model summary from trainHVT().

hvt.results$model_info$input_parameters

## $input_dataset
## [1] "160000 Rows & 3 Columns"
## 
## $n_cells
## [1] 100
## 
## $depth
## [1] 1
## 
## $quant.err
## [1] 0.1
## 
## $normalize
## [1] TRUE
## 
## $distance_metric
## [1] "L1_Norm"
## 
## $error_metric
## [1] "max"
## 
## $quant_method
## [1] "kmeans"
## 
## $diagnose
## [1] FALSE
## 
## $projection.scale
## [1] 10
## 
## $hvt_validation
## [1] FALSE
## 
## $train_validation_split_ratio
## [1] 0.8

Now, Let’s plot the Voronoi tessellation for 100 cells.

plotHVT(
  hvt.results,
  centroid.size = c(0.6),
  plot.type = '2Dhvt',
  cell_id = FALSE)

Figure 2: The Voronoi tessellation for layer 1 shown for the 100 cells in the dataset ’Lorenz attractor’

To understand how cell IDs are distributed across the map, we again plot Voronoi tessellation with cell_id = TRUE.

plotHVT(
  hvt.results,
  centroid.size = c(0.6),
  plot.type = '2Dhvt',
  cell_id = TRUE)

Figure 3: The Voronoi tessellation for layer 1 shown for the 100 cells in the dataset ’Lorenz attractor’ with Cell ID

6. Scoring

Now, once we have built the model, let us try to score using the model.

NOTE: we are using the training dataset plus the testing dataset i.e., the raw dataset here for the demo purpose.

set.seed(240)
scoring <- scoreHVT(dataset,
                    hvt.results,
                    child.level = 1)

Let’s see the output of scoreHVT() which has scored the cells for all the data points in the dataset. For the sake of brevity, we will only show the first 100 rows.

In the displayTable function, the ‘columnName’ argument takes any of the column names that have to be highlighted greater than the given ‘value’ based on the provided data.

displayTable(scoring$scoredPredictedData)

Let’s look at the scored model summary from scoreHVT()

The model info displays five attributes which are explained below:

1. Dimension of the dataset that is used in scoring.
2. The range of ‘Quantization Error’ after scoring of all the data points of the given dataset. (from minimum to maximum)
3. The value of Mean Absolute Deviation(MAD) Threshold.
4. The value of how many data points are anomalous. If the quantization error of a data point is greater than the ‘mad.threshold’, then it is considered as anomaly.
5. The value of how much cells have anomaly. If a cell has even 1 anomalous data point, the cell will be considered as anomaly_cell.

scoring$model_info$scored_model_summary

## $input_dataset
## [1] "200000 Rows & 5 Columns"
## 
## $scored_qe_range
## [1] "4e-04 to 0.3021"
## 
## $mad.threshold
## [1] 0.2
## 
## $no_of_anomaly_datapoints
## [1] 1353
## 
## $no_of_anomaly_cells
## [1] 17

Now, let’s merge back the U and t columns from the dataset to the scoring output and prepare the dataset for ‘Temporal Analysis’ Functions.

temporal_data <- cbind(scoring$scoredPredictedData, dataset[,c(4,5)]) %>% select(Cell.ID,t)

7. Timeseries plot with State Transitions

Let’s comprehend the function plotStateTransition which is used to create a time series plotly object.

plotStateTransition(
       df,
       sample_size,
       line_plot,
       cellid_column,
       time_column)

• df - A data frame contains Cell ID and Timestamps.

• sample_size - A numeric value to specify the sampling value which ranges between 0.1 to 1. The highest value 1, outputs a plot with the entire dataset. Sampling of data takes place from the last to the first.

• line_plot - A Logical value. If TRUE, the output will be a timeseries plot with a line connecting the states according to the sample_size. If FALSE, a timeseries plot but without a line based on the sample_size will be the output.

• cellid_column - A character string specifying the column name of the Cell ID.

• time_column - A character string specifying the column name of the time stamp.

plotStateTransition(df = temporal_data, 
                    cellid_column = "Cell.ID", 
                    time_column = "t",
                    sample_size = 0.2)

For the demo of ‘sample_size’ argument, we are replicating the same plot above with the entire dataset.

plotStateTransition(df = temporal_data, 
                    cellid_column = "Cell.ID", 
                    time_column = "t",
                    sample_size = 1)

8. Transition probability tables

getTransitionProbability(
        df, 
        cellid_column, 
        time_column)

• df - A data frame contains Cell ID and Timestamps.

• cellid_column - A character string specifying the column name of the Cell ID.

• time_column - A character string specifying the column name of the time stamp.

trans_table <- getTransitionProbability(df = temporal_data, 
                                        cellid_column = "Cell.ID", 
                                        time_column = "t")

NOTE: The output is stored as a nested list. For the purpose of demo, here we are displaying it as data frame with the first 10 rows.

combined_df <- do.call(rbind, trans_table)
displayTable(head(combined_df,10))

• Current_State: The cell (out of 100 cells given in model training) in which the datapoint resides at a given time (t).

• Next_State: The cell (out of 100 cells given in model training) to which the datapoint moves at the next time unit (t+1).

• Relative_Frequency: The number of times that the datapoint moves from that Current_State to that Next_State.

• Transition_Probability: The probability calculated from the Relative_Frequency. Individual Relative_Frequency divided by the total of Relative_Frequency for a particular Current_State.

9. Reconciling transition probability using markovchain package

Reconciling the transition probability of the current states to the next states manually and markovchain function considering self states and without self states.

reconcileTransitionProbability(
                df, 
                hmap_type = "All", 
                cellid_column, 
                time_column)

• df - A data frame contains Cell ID and Time stamps.

• hmap_type - A character string, if set to ‘without_self_state’, reconciliation plots for manual and markovchain for the highest transition probability excluding the self-state are given as output. If set to ‘self_state’, reconciliation plots for manual and markovchain for the highest transition probability considering the self-state is given as output and if set to ‘All’, plots including and excluding self-state are given as output.

• cellid_column - A character string specifying the column name of the Cell ID.

• time_column - A character string specifying the column name of the time stamp.

reconcile_plots <- reconcileTransitionProbability(df = temporal_data, 
                                                  hmap_type = "All", 
                                                  cellid_column = "Cell.ID",
                                                  time_column = "t")

Reconciliation plots of transition probability with self_state

The transition probability of one state staying in the same state is calculated using manual calculations and the markovchain function is plotted for comparison. The darker diagonal cells indicate higher probabilities of cells staying in the same state.

reconcile_plots[[1]]

Reconciliation table of transition probability with self-state

displayTable(reconcile_plots[[2]], limit = 217)

Reconciliation plots of transition probability without self-state

The transition probability of one state moving to the next state is calculated using manual calculations and the markovchain function and plotted for comparison. From all the next state transitions, the one with a higher probability is selected.

reconcile_plots[[3]]

Reconciliation table of transition probability without self-state

displayTable(reconcile_plots[[4]], limit = 217)

10. Animated Flowmaps

plotAnimatedFlowmap(
         hvt_model_output, 
         transition_probability_df, 
         df, 
         animation = "All", 
         flow_map = "All", 
         fps_state,
         fps_time,
         time_duration,
         state_duration,
         cellid_column, 
         time_column )

• hvt_model_output - The list object, which is the output from the trainHVT function.

• transition_probability_df - The probability list, which is the output from the getTransitionProbability function.

• df - A data frame contains Cell ID and Time stamps.

• animation - A character string, if set to ‘time_based’, an animation in which a red dot moves along cells according to the timestamp will be displayed. If set to ‘state_based’, an arrow animation based on the highest state excluding self-state will be displayed. If set to ‘All’, both the animation will be displayed. If set to NULL none will be displayed.

• flow_map - A character string. If set to ‘self_state’, the plot which shows the self-state by circles will be displayed. More the circle size, more probability of the cell stays in the same cell. If set to ‘without_self_state’, the plot which shows the next state by arrows will be displayed. The arrow head points the next cell to go from the cell which is pointed by the arrow tail. If set to ‘All’, two flow maps will be displayed. If set to NULL none will be displayed.

• fps_time - A numeric value indicating the frames per second for time transition animation. (Must be a numeric value and a factor of 100). Default value is 1.

• fps_state - A numeric value indicating the frames per second for state transition animation. (Must be numeric value and a factor of 100). Default value is 1.

• time_duration - A numeric value indicating the duration of the gif for time transition animation. Default value is 2.

• state_duration - A numeric value indicating the duration of the gif for state transition animation. Default value is 2.

• cellid_column - A character string specifying the column name of the Cell ID.

• time_column - A character string specifying the column name of the Time stamp.

flowmap_plots <- plotAnimatedFlowmap(hvt_model_output = hvt.results,
                                     transition_probability_df =trans_table,
                                     df = temporal_data, 
                                     animation = NULL , flow_map = 'All',
                                     fps_time = 30, fps_state =  5, 
                                     time_duration = 180,state_duration = 20,
                                     cellid_column = "Cell.ID", time_column = "t")
#> [1] "'animation' argument is NULL"

1. Flow map: Highest transition probability including self-state

The Circle size around the cell’s centroid represents self-state probability. More size, more probability of staying in the same cell.

flowmap_plots[[1]]

2. Flow map: Highest transition probability excluding self-states: Arrow size represents transition probability

The arrow size represents the Probability of the data to move to the next state. And the arrow directions point out to which cell it is moving next.

flowmap_plots[[2]]

3. Flow map animation: Highest state transition probabilities including self-state

The red point moves through the cells as per the time stamp and blinks mean that its stays in the same cell for a particular amount of time which can be referred from the sub-header in the gif.

#flowmap_plots[[3]] 

4. Flow map animation: Highest state transition probabilities excluding self-states

The arrow moves from the current cell to the next cell and its length refers to the transitioning probabilities of the next state without including self state.

#flowmap_plots[[4]]